#Plots produced with R, version 2.11
# Development Core Team (2010). R: A language and environment
#  for statistical computing. R Foundation for Statistical
 # Computing, Vienna, Austria. ISBN 3-900051-07-0, URL
  #http://www.R-project.org.


rm(list=ls(all=TRUE))
library(foreign)
library(lattice)

#Set working directory to directory holding data files 
data <- read.dta("replication_5-5-10.dta")

##############
###FIGURES 3 & 4###
##############
data.fs <- data[,c(1:2,14,19,25)]
data.fs <- data.fs[data$frac_rural>.1,]
data.fs <- na.omit(data.fs)
data.fs$rain_monthly <- data.fs$rain_monthly - ave(data.fs$rain_monthly, data.fs$code)
data.fs$rain_monthly <- data.fs$rain_monthly - ave(data.fs$rain_monthly, data.fs$year)
data.fs$rain_monthly_noabs <- data.fs$rain_monthly_noabs - ave(data.fs$rain_monthly_noabs, data.fs$code)
data.fs$rain_monthly_noabs <- data.fs$rain_monthly_noabs - ave(data.fs$rain_monthly_noabs, data.fs$year)
data.fs$ag_income <- data.fs$ag_income - ave(data.fs$ag_income, data.fs$code)
data.fs$ag_income <- data.fs$ag_income - ave(data.fs$ag_income, data.fs$year)

#function for loess model fitting and bootstrapping
smooth.bt <- function(xy, bw, nloc=50, nbt=100,lwd=3,col="blue",
                          smooth.func=function(xy,bw) lowess(xy,f=bw)) {
      x.boot <- matrix(ncol=nbt,nrow=dim(xy)[1])
      y.boot <- matrix(ncol=nbt,nrow=dim(xy)[1])
      for(i in 1:nbt) {
        sel <- sample(1:dim(xy)[1], replace=T)
        fit <- smooth.func(xy[sel,],bw=bw)
        x.boot[,i] <- fit$x
        y.boot[,i] <- fit$y
        #lines(fit, col="gray50")
        print(i)
      }
      bs.pts <- smooth.func(xy, bw=bw)
      bs.x.pts <- bs.pts$x
      bs.se <- apply(y.boot,1,quantile,probs=c(.05,.95))
      list(bs.pts=bs.pts, bs.se = bs.se)
          }

lowess.first <- smooth.bt(xy=cbind(data.fs$rain_monthly_noabs, data.fs$ag_income),nbt=50,bw=.6)

pdf("firststage.pdf")
plot(c(min(data.fs$rain_monthly_noabs),max(data.fs$rain_monthly_noabs)),c(-.15,.045),type = "n", xlab = "Rainfall",ylab = "Agricultural Income")
lines(lowess.first$bs.pts,lwd=2,col=1)
lines(lowess.first$bs.pts$x,lowess.first$bs.se[1,],lty=2,col="black")
lines(lowess.first$bs.pts$x,lowess.first$bs.se[2,],lty=2,col="black")
text(-2.5,0.045,"Bandwidth=0.6")
dev.off()   


lowess.first.abs <- smooth.bt(xy=cbind(data.fs$rain_monthly, data.fs$ag_income),nbt=50,bw=.6)

pdf("firststage_abs.pdf")
plot(c(min(data.fs$rain_monthly),max(data.fs$rain_monthly)),c(-.1,.09),type = "n", xlab = "Rainfall",ylab = "Agricultural Income")
lines(lowess.first.abs$bs.pts,lwd=2,col=1)
lines(lowess.first.abs$bs.pts$x,lowess.first.abs$bs.se[1,],lty=2,col="black")
lines(lowess.first.abs$bs.pts$x,lowess.first.abs$bs.se[2,],lty=2,col="black")
text(-.3,0.09,"Bandwidth=0.6")
dev.off()   





##############
###FIGURE 4###
##############
data.3d <- data[,c(1:2,6,25,12)]
#remove missing data
data.3d <- na.omit(data.3d)

#Take out fixed effect
data.3d$occs <- data.3d$occs - ave(data.3d$occs, data.3d$code)
data.3d$occs <- data.3d$occs - ave(data.3d$occs, data.3d$year)
data.3d$rain_monthly <- data.3d$rain_monthly - ave(data.3d$rain_monthly, data.3d$code)
data.3d$rain_monthly <- data.3d$rain_monthly - ave(data.3d$rain_monthly, data.3d$year)

#estimate loess model
loess.3d <- loess(occs~land_gini*rain_monthly,data=data.3d,span=.6,degree=1)

#Create dataset to hold predicted values
plot.3d.data <- data.frame(rain_monthly = seq(min(data.3d$rain_monthly),max(data.3d$rain_monthly),len=25))
plot.3d.data$land_gini <- seq(min(data.3d$land_gini),max(data.3d$land_gini),len=25)
plot.3d.data <- expand.grid(rain_monthly=plot.3d.data$rain_monthly, land_gini=plot.3d.data$land_gini)
#Create predicted number of occupations
plot.3d.data$occs <- matrix(predict(loess.3d, newdata = plot.3d.data), 625, 1)

#Create plot
trellis.par.set("axis.line",list(col="transparent"))
plot.3d <- wireframe(occs~rain_monthly*land_gini,
                data=plot.3d.data,
                scales=list(arrows=FALSE,
                  cex=.8,col="black",font=3,tck=1),
                drape=FALSE,
                colorkey=FALSE,
                screen=list(z=-305,x=-65),
                xlab= list(label="Rainfall", cex=.75),
                ylab= list(label="Land Gini",cex=.75),
                zlab= list(label="Invasions",cex=.75),
                pretty=TRUE
                )

pdf("invasions_rain_gini.pdf")
plot.3d
dev.off()   


##############
###FIGURE 5###
##############

#Load cross-section dataset
data_crossection <- read.dta("replication_cs_5-5-10.dta")
data_crossection <- data_crossection[data_crossection$frac_rural>.1,]
data_crossection <- data_crossection[is.na(data_crossection$landgini)==FALSE & is.na(data_crossection$occs)==FALSE,]

#take out micro-region fixed effects
data_crossection$FixOccs <- data_crossection$occs - ave(data_crossection$occs,data_crossection$microregion)
data_crossection$FixGini <- data_crossection$landgini - ave(data_crossection$landgini,data_crossection$microregion)

#Fit loess model
loess.cross <- loess(data_crossection$FixOccs~data_crossection$FixGini, span =.75)
#Bootstrap
loess.boot <- matrix(nrow=nrow(data_crossection),ncol=100)
for(i in 1:100){
  boot.index <- sample(1:4062, 4062, replace=TRUE)
  loess.boot[,i] <- loess(data_crossection$FixOccs[boot.index]~data_crossection$FixGini[boot.index], span =.75)$fitted[order(data_crossection$FixGini[boot.index])]
}
loess.ci <- apply(loess.boot, 1, quantile,probs=c(.05,.95))

#Plot
pdf("crossSectionLoess.pdf")
plot(data_crossection$FixGini,data_crossection$FixOccs,cex=.001,type="n",ylim=c(-2.7,3.8),xlab="Land Gini",ylab="Occupations")
lines(spline(data_crossection$FixGini,loess.cross$fitted),lwd=2,col="darkgreen")
lines(sort(data_crossection$FixGini),loess.ci[1,],lty=2,lwd=2,col="darkgreen")
lines(sort(data_crossection$FixGini),loess.ci[2,],lty=2,lwd=2,col="darkgreen")
text(-.33,3.8,"Bandwidth=0.75")
dev.off()   
